Sesame, Door Opens.
As well-known history, Sir Newton used the tools of Calculus to study mechanics.
Let \begin{eqnarray} T_0 & : & \text{the initial temperature of a body at } t_0\ T (t) & : & \text{the temperature of the body at time $t$,}\ T_m & : & \text{the temperature of certain a medium and $T_m < T_0$.} \end{eqnarray} Now put the body into the medium. The cooling law states that the rate which the body cools at any given instant is proportional to the difference between its at that instant and the the temperature $T_m$ of the surrounding medium, i.e. $$ \frac{T(t+\Delta t)-T(t)}{\Delta t}\sim T(t)-T_m $$
\begin{eqnarray} & T' (t) & = k (T (t) - T_m)\ & T (t_0) & = T_0\ \Longrightarrow & \int \frac{d T}{T - T_m} & = \int k d t\ \Longrightarrow & \ln (T - T_m) & = k t + C\ \Longrightarrow & T & = T_m + c e^{k t}\ \Longrightarrow & T & = T_m + (T_0 - T_m) e^{ k (t-t_0)} \end{eqnarray} by the initial conditions.
from sympy import *
from sympy.abc import t,x,k,N
init_printing()
var('C1 C2 C3')
y=Function("y")
def Dsol(func,t0,initial):
"""
ODE solver
"""
y=Function("y")
yp=Derivative(y(t),t)
sol=dsolve(yp-func,y(t))
cc=solve(sol.rhs.subs(t,t0)-initial,C1)
return sol.rhs.subs(C1,cc[0])
var("L t0 T0 Tm")
func=k*(y(t)-Tm)
Dsol(func,t0,T0)
A dead body was found around Carven Garden, London. The Holmes and Dr. Watson arrived the scene at 9 : 00 A.M. and immediately recorded the temperature of the body as $30.1^\circ$C. One hour later they noted that the body's temperature had dropped to $29.2^\circ$C, and during the hour, room temperature had remained at $20^\circ$ˆ˜C. Assume that the body's temperature obeys Newton's law of cooling, determine when death occurred. (Normal body temperature is $37^\circ$C.)
Which is the following is right? a) AM 2hr and 09min b) AM 3hr and 26min c) AM 4hr and 47min
The path an object takes when chasing after another object in the most effective way.
In Francis Godwin’s story, The Man in the Moone, an astronaut harnesses a wedge of 25 swans and flies to the moon. The swans fly at a constant rate and always head toward the moon, in accordance with their annual migration. Thus, the trajectory from the earth to the moon is not a straight line, but is a pursuit curve. In the story, Godwin says that the flight to the moon takes twelve days, whereas the return trip, which follows a straight line, takes eight days days.
References:

0. Basics of Numerical: from Continuous to discrete
1. Basics of Physics;
2. Pictures describing the tractories of observerved objects;
3. Animation of whole the processes;
4. Simple Interaction by passing the Parameter to System.
1. Python >3.x
2. numpy: data handling
3. matplotlib-1.3.1: making pictures and animation,
4. IPython, JSAnimation, IPyWidget: make results well display and manipulation from HTML-based browsers.
Consider the following initial value problem:
\begin{array}{cccl}
&x'(t)&=&x(t), \cr
&x(0)&=&1 \cr
\Rightarrow &x &=& \exp(t).
\end{array}
The respective numerical apprroximation, $$x'(t) \sim \frac {x(t+\vartriangle t)-x(t)}{\vartriangle t}.$$
Now partition $[0,t]$ into $n$ equally-like nonoverlapping sub-intervals as follows:
\begin{array}{cccccc} x: & x_0 = 1 & x_1 & x_2 & \cdots & x_{n-1} & \quad x_n \cr & \bot \ldots & ... \perp ... & ... \perp ... & \ldots \ldots . & ... \perp ... & \ldots \bot\cr t: & 0\quad & \vartriangle t & 2\vartriangle t & \cdots & (n-1)\vartriangle t & \quad n\vartriangle t=t \end{array}where $x_i=x(i\vartriangle t),$ and $\vartriangle t=\frac{t}{n}$.
This induces the following difference scheme: \begin{array}{cccl} &x'&=& x \cr &\frac{x_{i+1}-x_i}{\vartriangle t} &\sim& xi \cr \Rightarrow &x{i+1} &\sim& (1+\vartriangle t)xi &\sim& (1+\vartriangle t)^2x{i-1} \ &&\sim& (1+\vartriangle t)^{i+1} x_0 \end{array} for $i=1,2,\cdots$.
$a^\circ$) Estimate the value of Euler number,$e\sim2.71828\cdots$
$e=\exp(1)$ could be evaluated by
\begin{array}{ccl}
x_0 &=& e(0), \cr
x_1 &=& e(\vartriangle t) \cr
&\sim& (1+\vartriangle t)x_0, \cr
x_2 &=& e(2\vartriangle t) \cr
&\sim& (1+\vartriangle t)x_1 \cr
&\sim& (1+\vartriangle t)^2 x_0 \cr
&\vdots \cr
x_n &\sim& (1+\vartriangle t)^n x_0
\end{array}
, which refers to the following partition for $\exp(x)$ from 0 to 1:
\begin{array}{cccccc}
x: & x_0 = 1 & x_1 & x2 & \cdots & x{n-1} & x_n = e\cr
&\bot \ldots & ... \perp ... & ... \perp ... & \ldots \ldots . & ... \perp
... & \ldots \bot\cr
t: & 0\quad & \frac{1}{n} & \frac{2}{n} & \cdots & \frac{n-1}{n} & \quad1
\end{array}
where $\vartriangle t=\frac{1}{n}$. And the approximation for Euler number, $e$, could be estimeted by the following: $$ e=x_n\sim \left(1+\frac{1}{n}\right)x_{n-1}\sim \left(1+\frac{1}{n}\right)^2x_{n-2}\sim\dots \sim \left(1+\frac{1}{n}\right)^n e(0)=\left(1+\frac{1}{n}\right)^n$$
n=100
a=0;b=1;
dt=(b-a)/float(n)
e_value=(1+dt)**n
print 'the approximation of e is ', e_value
from IPython.display import clear_output
import sys
import time
def sim(n):
a=0;b=1;
dt=(b-a)/float(n)
e_value=(1+dt)**n
print "the approximation, ( 1 + 1/",n,")^",n,", of e is " ,e_value
for i in range(1,11):
time.sleep(0.5)
clear_output()
n=10*i
sim(n)
sys.stdout.flush()
%matplotlib inline
import numpy as np
from numpy import exp
import matplotlib.pylab as plt
"""
from pylab import *
rcParams['figure.figsize'] = 4,4
"""
n=11
a=0.;b=1.
dx=(b-a)/float(n-1)
x=np.linspace(a,b,n)
y=np.ones(n)
for i in np.arange(n):
y[1:]=(1+dx)*y[:-1]
P2=plt.plot(exp(x),'b-',y,'r*')
$b^\circ$) Logistic Model (or called gassip model, learning model,...)
\begin{align} y'&=&a y (M-y) \cr y(0) &=& y_0 \end{align}
where $a,M, y_0$ are constants and $M$ is the maximal resources that environment could afford.
The scheme runs as follows: \begin{array}{ccl} &\frac{y_{n+1}-y_n}{\vartriangle t} &=& a y_n (M-yn) \cr \Rightarrow &y{n+1} &=& y_n +a y_n (M-y_n)\vartriangle t \end{array}
The solution curve is like $S$-shape:
a=1/5.
M=20
t0=0;t1=1
n=30
u0=5
t=np.linspace(t0,t0+t1,n+1)
dt=(t1-t0)/float(n)
u=np.zeros(np.size(t))
u=u0*np.ones(np.size(u))
u_h=M/2.*np.ones(np.size(u))
for i in np.arange(len(u)-1):
"""
i.e.
u[i+1] = u[i]+a*u[i]*(M-u[i])*dt
"""
u[1:]=u[:-1]+a*u[:-1]*(M-u[:-1])*dt
plt.plot(t,u,t,u_h,'r--')
Levels=[1,2,5,10,13,15,17,25,22,21]
P=[]
for level in Levels:
u=np.zeros(np.size(t))
u[0]=level
for i in np.arange(len(u)-1):
u[1:]=u[:-1]+a*u[:-1]*(M-u[:-1])*dt
P+=plt.plot(t,u)
P+=plt.plot(u_h,'r--')
P+=plt.plot(2*u_h,'b--')
plt.axis([0,1,0,26])
#show(P)
Newton's law of cooling states that the rate change of temperate of heated object is propotional to the difference between its own temperature and the ambient temperature. $$ \mathbf{y(t)=k(y(t)-y_m)}$$ where $y(0)=y_0$ is initial temperature of temperature, $y'(t)$ is the temperature at time $t$, and $y_m$ is ambient temperature. Suppose that a corpse was discovered in a room and its temperature was 32°C. The temperature of the room is kept constant at 22°C. Three hours later the temperature of the corpse dropped to 27°C.
Example Why is there no people died of being hit by falling rain?
According Newton's Law, the rain drop will fall down as a free body and be formulated by the following: $$\ddot{\mathbf{x}} = \mathbf{g}-r|\dot{\mathbf{x}}|^2 $$
import numpy as np
import scipy.integrate as sp
import matplotlib.pyplot as plt
%matplotlib inline
# The initial position and velocity is (0,0) (or * above sea level)
x0 = np.zeros(2)
g=9.81
r=1.
def f(x,t, r):
# X has twocomponents: X=[x, x'].
x, xdot = x[:1], x[1:]
# We compute the second derivative x'' of x.
xdotdot = g-r * xdot *xdot
# We return X'=[x', x''].
return np.r_[xdot, xdotdot]
t = np.linspace(0., 10., 100)
x = sp.odeint(f, x0, t, args=(r,))
from numpy import sqrt
plt.plot(t,sqrt(g)*np.ones(np.size(t)))
plt.plot(t[::2],x[::2,1], 'o-', mew=1, ms=8, mec='w',mfc="r")
from sympy import Symbol , dsolve , Function , Derivative ,symbols, Symbol,solve, simplify,init_printing,integrate,exp
#from sympy.abc import a,r,t
init_printing()
t = Symbol("t",positive=True)
r,x = symbols("r x",positive=True)
v=Function('v')
g=9.8
#r=1
f_=Derivative(v(t),t)-g+r*v(t);print f_
sol=dsolve(f_,v(t));print sol.rhs
Suppose that a dog, moving at constant speed while always pointing towards its master.
Consider:
The motion could be represented simply by the vector form:
$$\vec{v_D}=|\vec{v_D}| \cdot \frac{\vec{M}-\vec{D}}{||\vec M-\vec D||}$$
Then the coordinates of position is as follows:
%matplotlib inline
import numpy as np
from numpy import pi,sin,cos,size
from numpy.linalg import norm
import matplotlib.pyplot as plt
from matplotlib import animation
from JSAnimation import IPython_display
from ipywidgets import interact, IntSlider#,Radio
from IPython.display import display, clear_output
Dpos=np.array([[0,0]])
Mpos=np.array([[0,10]])
dt=0.1
Mvel=np.array([5,0])
#while i<40:
#while norm(Dpos[-1]-Mpos[-1])>0.01:
while Dpos[-1][1]<Mpos[0][1]:
Dvel=10*(Mpos[-1]-Dpos[-1])/norm(Mpos[-1]-Dpos[-1])
Mpos1 = Mpos[-1]+Mvel*dt
Dpos2 = Dpos[-1]+Dvel*dt
Mpos=np.vstack([Mpos,Mpos1])
Dpos=np.vstack([Dpos,Dpos2])
#i=i+1
figM = plt.figure(figsize=(6,6))
ax = figM.add_subplot(111)
ax.set_title("Dog Chasing",fontsize=14)
plt.xlim([0,12])
plt.ylim([0,12])
ax.scatter(Mpos[:,0],Mpos[:,1], c='r' )
ax.scatter(Dpos[:,0],Dpos[:,1],c='g')
def DogTrack(Vel):
#while i<40:
#while norm(Dpos[-1]-Mpos[-1])>0.01:
Dpos=np.array([[0,0]])
Mpos=np.array([[0,10]])
i=1
#while Dpos[-1][1]<Mpos[0][1] or i<40:
while i<40:
Dvel=Vel*(Mpos[-1]-Dpos[-1])/norm(Mpos[-1]-Dpos[-1])
Mpos1 = Mpos[-1]+Mvel*dt
Dpos2 = Dpos[-1]+Dvel*dt
Mpos=np.vstack([Mpos,Mpos1])
Dpos=np.vstack([Dpos,Dpos2])
i=i+1
figM = plt.figure(figsize=(9,7))
#ax = figM.add_subplot(211)
plt.title("Dog Chasing",size=14)
plt.xlim([0,18])
plt.ylim([0,14])
plt.scatter(Mpos[:,0],Mpos[:,1], c='r' )
plt.scatter(Dpos[:,0],Dpos[:,1],c='g')
plt.text(1, 11.5,
"Master Velocity={0}\nDog Velocity={1:.2f}".format(5,Vel),
size=14, color='gray')
#return figM
interact(DogTrack,Vel=IntSlider(min=2, max=10, step=1))
Four cats are in the four corners a squareroom of unitary base. At time t=0, they are starting to chase the tails of others. Same as the dog tracking problem:
The case of the cat at the left lower corner runs twice faster than others is simulated as follows:
$$vel_i=\frac{\vec{P}_i-\vec{P}_{i-1}}{||\vec{P}_i-\vec{P}_{i-1}||}$$$$ \vec P_{i}(t_{n+1})=\vec P_{i}(t_n)+vel_{i}(t_n)\times d t$$Cpos1=np.array([[0,0]])
Cpos2=np.array([[1,0]])
Cpos3=np.array([[1,1]])
Cpos4=np.array([[0,1]])
dt=0.1
i=1
#while i<40:
while norm(Cpos1[-1]-Cpos2[-1])>0.01:
v1=2*(Cpos2[-1]-Cpos1[-1])
v2=2*(Cpos3[-1]-Cpos2[-1])
v3=2*(Cpos4[-1]-Cpos3[-1])
v4=3*(Cpos1[-1]-Cpos4[-1])
Cpos11 = Cpos1[-1]+v1*dt
Cpos21 = Cpos2[-1]+v2*dt
Cpos31 = Cpos3[-1]+v3*dt
Cpos41 = Cpos4[-1]+v4*dt
Cpos1=np.vstack([Cpos1,Cpos11])
Cpos2=np.vstack([Cpos2,Cpos21])
Cpos3=np.vstack([Cpos3,Cpos31])
Cpos4=np.vstack([Cpos4,Cpos41])
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.set_title("Four Cats Pursuit Curve",fontsize=14)
plt.xlim([0,1])
plt.ylim([0,1])
ax.scatter(Cpos1[:,1],Cpos1[:,0], c='y' )
ax.scatter(Cpos2[:,1],Cpos2[:,0],c='g')
ax.scatter(Cpos3[:,1],Cpos3[:,0],c='b')
ax.scatter(Cpos4[:,1],Cpos4[:,0],c='r')
anim = plt.figure(figsize=(6,6))
ax = anim.add_subplot(111)
ax.set_title("Four Cats Pursuit Curve",fontsize=14)
plt.xlim([0,1])
plt.ylim([0,1])
def init():
return ax.scatter(Cpos1[0,1],Cpos1[0,0], c='y' ), \
ax.scatter(Cpos2[0,1],Cpos2[0,0],c='g'), \
ax.scatter(Cpos3[0,1],Cpos3[0,0],c='b'), \
ax.scatter(Cpos4[0,1],Cpos4[0,0],c='r')
def animate(i):
return ax.scatter(Cpos1[0:i,1],Cpos1[0:i,0], c='y' ), \
ax.scatter(Cpos2[0:i,1],Cpos2[0:i,0],c='g'), \
ax.scatter(Cpos3[0:i,1],Cpos3[0:i,0],c='b'), \
ax.scatter(Cpos4[0:i,1],Cpos4[0:i,0],c='r')
animation.FuncAnimation(anim, animate, init_func=init, frames=size(Cpos1[:,1]), interval=size(Cpos1[:,1]))
Suppose that one is runs faster than others, Then ... This is nothing to be worry because of the trusted simulation scheme doing well.
Cpos1=np.array([[0,0]])
Cpos2=np.array([[1,0]])
Cpos3=np.array([[1,1]])
Cpos4=np.array([[0,1]])
#while i<40:
while norm(Cpos1[-1]-Cpos2[-1])>0.01:
v1=(Cpos2[-1]-Cpos1[-1])
v2=(Cpos3[-1]-Cpos2[-1])
v3=(Cpos4[-1]-Cpos3[-1])
v4=(Cpos1[-1]-Cpos4[-1])
Cpos11 = Cpos1[-1]+3.*v1*dt
Cpos21 = Cpos2[-1]+v2*dt
Cpos31 = Cpos3[-1]+v3*dt
Cpos41 = Cpos4[-1]+v4*dt
Cpos1=np.vstack([Cpos1,Cpos11])
Cpos2=np.vstack([Cpos2,Cpos21])
Cpos3=np.vstack([Cpos3,Cpos31])
Cpos4=np.vstack([Cpos4,Cpos41])
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.set_title("Four Cats Pursuit Curve",fontsize=14)
plt.xlim([0,1])
plt.ylim([0,1])
ax.scatter(Cpos1[:,1],Cpos1[:,0], c='y' )
ax.scatter(Cpos2[:,1],Cpos2[:,0],c='g')
ax.scatter(Cpos3[:,1],Cpos3[:,0],c='b')
ax.scatter(Cpos4[:,1],Cpos4[:,0],c='r')
Suppose that there are three cats stand at the corners of the equilateral triangle. Guess the trajectories about the cats chansing others' tails.
If they run wih the different speed with respect to others, what would happen?
In 1968, Apollo 8, the first manned craft to circle the moon, reached the moon in 3.5 days, and returned in 2.5 days.
How fast does the Moon travel around Earth?
The Moon orbits Earth at a speed of 2,288 miles per hour (3,683 kilometers per hour). During this time it travels a distance of 1,423,000 miles (2,290,000 kilometers).
$$ 1 \text{ Lunar distance }=1 LU\sim 384,000 \text{ km }\\ p \text{ period } \sim 29.5 \text{ days}\\ R_{\text{earth}}~0.067 LU, r_{\text{moon}}~0.00453 LU \\ M(t): \text{ position of Moon }= (M_1(t),M_2(t))=\left(\cos\left(\frac{2\pi t}{p}\right), \sin\left(\frac{2\pi t}{p}\right)\right)\\ (x,y): \text{ position of spacecraft} $$where $(x,y)$ is coordinate of spacecraft and $(M_1,M_2)$ is the position of Lunar.
Mpos=np.array([[1,0]])
Epos=np.array([[0,0]])
p=29.5
c=0.15
dt=0.1
i=1
while i<400:
#while norm(Dpos[-1]-Mpos[-1])>0.01:
#while Dpos[-1][1]<Mpos[0][1]:
Mpos1=np.array([cos(2*pi*i*dt/p),sin(2*pi*i*dt/p)])
Evel=c*(Mpos[-1]-Epos[-1])/norm(Mpos[-1]-Epos[-1])
#print norm(Mpos[-1]-Epos[-1])
Epos2 = Epos[-1]+Evel*dt
#print Evel,Epos2
Mpos=np.vstack([Mpos,Mpos1])
Epos=np.vstack([Epos,Epos2])
i=i+1
figM = plt.figure(figsize=(6,6))
ax = figM.add_subplot(111)
ax.set_title("Moon Landing",fontsize=14)
plt.xlim([-1.1,1.1])
plt.ylim([-1.1,1.1])
ax.scatter(Mpos[:,0],Mpos[:,1], c='r' )
ax.scatter(Mpos[-1,0],Mpos[-1,1], c='b',s=200)
ax.scatter(Epos[-1,0],Epos[-1,1], c='b',s=80)
ax.plot([Mpos[-1,0],Epos[-1,0]],[Mpos[-1,1],Epos[-1,1]],'b-->')
#ax.quiver(Mpos[-1,0]-0.2, Mpos[-1,1]-0.2, (Mpos[-1,0]-Epos[-1,0]), (Mpos[-1,1]-Epos[-1,1]), pivot='middle',\
# headwidth=1, headlength=1,color="red")
ax.scatter(Epos[:,0],Epos[:,1],c='g')
Epos[-1,0], Epos[-1,1], Mpos[-1,0]-Epos[-1,0], Mpos[-1,1]-Epos[-1,1]
anim = plt.figure(figsize=(6,6))
ax = anim.add_subplot(111)
ax.set_title("Moon Landing",fontsize=14)
ax.text(-0.45, -0.5,
"Velocity of Spacecraft={0:.2f}".format(c),
fontsize=14, color='gray')
plt.xlim([-1.1,1.1])
plt.ylim([-1.1,1.1])
def init():
return ax.scatter(Mpos[0,0],Mpos[0,1], c='r' ), \
ax.scatter(Epos[0,0],Epos[0,1],c='g')
def animate(i):
return ax.scatter(Mpos[0:4*i:4,0],Mpos[0:4*i:4,1], c='r' ), \
ax.scatter(Epos[0:4*i:4,0],Epos[0:4*i:4,1],c='g'), \
ax.plot([Mpos[4*i,0],Epos[4*i,0]],[Mpos[4*i,1],Epos[4*i,1]],'b--')
animation.FuncAnimation(anim, animate, init_func=init, frames=size(Mpos[::4,1]), interval=size(Mpos[::4,1]))
def MoonTraj(vec):
c=vec
Mpos=np.array([[1,0]])
Epos=np.array([[0,0]])
i=1
#while Dpos[-1][1]<Mpos[0][1] or i<40:
while i<800:
Mpos1=np.array([cos(2*pi*i*dt/p),sin(2*pi*i*dt/p)])
Evel=c*(Mpos[-1]-Epos[-1])/norm(Mpos[-1]-Epos[-1])
Epos2 = Epos[-1]+Evel*dt
Mpos=np.vstack([Mpos,Mpos1])
Epos=np.vstack([Epos,Epos2])
i=i+1
figM = plt.figure(figsize=(6,6))
#ax = figM.add_subplot(111)
plt.title("Moon Landing",size=14)
plt.xlim([-1.1,1.1])
plt.ylim([-1.1,1.1])
plt.scatter(Mpos[:,0],Mpos[:,1], c='r' )
plt.scatter(Epos[::4,0],Epos[::4,1],c='g')
plt.text(-0.5, -0.5,
"Velocity of Spacecraft={0:.2f}".format(c),
size=14, color='gray')
#return figM
from ipywidgets import FloatSlider
interact(MoonTraj,vec=FloatSlider(min=0.05, max=0.30, step=0.01))
If the velocity of projectile is not fast, it should fall down. Find out the minimal speed necessary to the moon.
Consider 2D ODE's as follows:
\begin{eqnarray} \frac{ d x}{ d t} &=& 2 x (1-x/2 ) - x y\\ \frac{ d y}{ d t} &=& 3 y (1 -y/3) - 2 x y \end{eqnarray}def f(x,y): return 2.*x*(1-x/2.) - x*y
def g(x,y): return 3.*y*(1.-y/3.) - 2.*x*y
It could describe how the system revolve:
fig1 = plt.figure(figsize=(9,6))
ax = fig1.add_subplot(1,1,1)
xmin=-0.5; xmax=3.5; ymin=-0.5; ymax=4.5
xmesh,ymesh=np.linspace(xmin,xmax,25),np.linspace(ymin,ymax,21)
plotwidth=5.0 # inches, horizontal distance between xmin and xmax on display
plotheight=3.9 # inches, vertical distance between ymin and ymax on display
alpha= plotwidth/(xmax-xmin); beta = plotheight/(ymax-ymin) # scaling factors
for xp in xmesh:
for yp in ymesh:
vx=f(xp,yp); vy=g(xp,yp) # direction vector v=(vx,vy)
ax.quiver(xp, yp, vx, vy, pivot='middle', headwidth=3, headlength=4,color="red")
plt.scatter((0,1,2,0),(3,1,0,0),c='b',marker='o',alpha=0.6)
ax.set_xlabel('x'); ax.set_ylabel('y')
ax.set_aspect(beta/alpha)
ax.set_xlim(xmin,xmax); ax.set_ylim(ymin,ymax)
plt.grid() # turn on a grid of light, dotted lines
plt.autoscale(enable=True,tight=True) # remove extra whitespace around slopefield
plt.title("Direction Field for\n $x'= 2x(1-x/2)-xy, y'= 3y(1-y/3)-2xy$", size=16)
There are four stabilities in system;
!jupyter nbconvert DE.ipynb